#include "../../headers/structsystem.h"
using namespace std;

inline void printStdArray6d(const StdArray6d* arr)
{
    for (int i = 0; i < 6; i++)
    {
        cout << (*arr)[i] << ",";
    }
    cout << endl;
}


int main()
{
    // create system
    StructSystem system = StructSystem(1);
    system.setJobname("Example1");
    double zeta = 0.0;
    cout << "Please enter a damping coefficient(>=0): " << endl;
    cin >> zeta;

    // solve parameters
    double endTime = 10.0;  // total time, s
    double h = 1.0e-3;      // time step, s
    double v0 = 10;         // initial velocity, m/s
    double mass = 10;       // mass, kg

    // create particles
    Particle* p1 = new Particle(1, 0, 0);
    system.addParticle(p1);
    // set mass properties
    StdArray6d m = {mass, 2*mass, 3*mass, 4*mass, 5*mass, 6*mass};
    p1->setLumpedMass(m);
    Eigen::Matrix3d Im;
    Im << 10, 1, 1,
         1, 10, 1,
         1, 1, 10;
    p1->setInertiaMass(Im);

    // activate dof of particle
    p1->activateDof("Ux");
    p1->activateDof("Uy");
    p1->activateDof("Uz");
    p1->activateDof("Rotx");
    p1->activateDof("Roty");
    p1->activateDof("Rotz");

    string path = system.workdir() + "/" + system.jobname() + "/model";
    system.saveModel(path);

    StdArray6d f = {10, 10, 10, 10, 10, 10};
    int nStep = ceil(endTime / h);
    int interval = ceil(0.1/h);
    for (int i = 0; i <= nStep; i++)
    {
        if (i == 0)
        {
            p1->setForce(f);
            system.solve(h, zeta, true);
        }
        else
        {
            system.solve(h, zeta);
            system.clearParticleForce();
            p1->setForce(f);
        }

        // save results
        if (i % interval == 0)
        {
            string path = system.workdir() + "/" + system.jobname() + "/" +
                          to_string(i*h);
            system.saveParticleResult(path);
        }
    }

    system.releaseContainers();
    return 0;
}
